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At the end of inflation, dynamical instability can rapidly deposit the energy of homogeneous 
cold inflaton into excitations of other fields. This process, known as preheating, is rather violent, 
inhomogeneous and non- linear, and has to be studied numerically. This paper presents a new code 
for simulating scalar field dynamics in expanding universe written for that purpose. Compared to 
available alternatives, it significantly improves both the speed and the accuracy of calculations, and 
is fully instrumented for 3D visualization. We reproduce previously published results on preheating 
in simple chaotic inflation models, and further investigate non-linear dynamics of the infiaton decay. 
Surprisingly, we flnd that the fields do not want to thermalize quite the way one would think. Instead 
of directly reaching equilibrium, the evolution appears to be stuck in a rather simple but quite 
inhomogeneous state. In particular, one-point distribution function of total energy density appears 
to be universal among various two-field preheating models, and is exceedingly well described by a 
lognormal distribution. It is tempting to attribute this state to scalar field turbulence. 

PACS numbers: 98.80.Cq, 05.10.-a 



I. INTRODUCTION 

The idea of inflation is a cornerstone of the modern 
theory of the early universe. According to inflationary 
paradigm, universe at early times undergoes a period of 
rapid (quasi-exponential) expansion, which wipes the ini- 
tial state of the universe clean while seeding the primor- 
dial inhomogeneities with quantum fluctuations gener- 
ated during expansion [l], Q . While universe is inflating, 
all of its energy sits in the homogeneous scalar field or 
condensate (known as inflaton), which is in a vacuum- 
like state with little entropy or particle excitations. But 
eventually inflation ends, and this energy has to be de- 
posited into excitations of other matter fields, starting 
the thermal history of the universe with a hot big bang. 

Decay of the infiaton can be very efficient if the fields 
experience dynamical instability at the end of infia- 
tion; such a stage became known as preheating. In 
most chaotic infiation models, oscillations of the infia- 
ton field can cause instability via parametric resonance 
S 01 H H ■ Although linear development of this in- 
stability can be understood analytically d, Q , it might 
be chaotic [l^l, and one needs to resort to numerical 
simulations to investigate non-linear dynamics that soon 
takes over [ll|, [l^, [HI, In hybrid inflation models 

p^ . in addition to parametric resonance p^ . one also 
has a tachyonic instability associated with the symmetry 
breaking [l7l|. dynamics of which has been explored in 

Non-equilibrium dynamics of preheating can lead to a 
multitude of interesting phenomena. Some of the topics 
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discussed in the literature are formation of topological 
defects [2ll . [22| , production of various particles (with ap- 
plications to baryo- and leptogenesis) [l^, 
possibility of primordial black hole formation 
generation of primordial magnetic fields [1^, |30|, and 
production of stochastic gravitational wave background 
[3l|, m, m, S, m, [3^ . Due to difliculties of dealing with 
non-linear evolution equations, most of these studies rely 
on numerical simulations. 

This paper describes a new code for simulating non- 
linear scalar field dynamics in expanding universe devel- 
oped to study preheating, called DEFROST, and the first 
results obtained with it. There are other codes avail- 
able for this purpose, most notably LATTICEEASY by 
Gary Felder and Igor Tkachev js^], and its parallel ver- 
sion CLUSTEREASY [H. Through the use of more 
advanced algorithms and careful optimization, the new 
code significantly improves both accuracy and perfor- 
mance achievable in simulations of preheating. An im- 
portant design goal has been the ease of visualization 
and analysis of the results, which is all too important for 
understanding dynamics of complex systems. 

We present results on preheating in a simple two- 
field chaotic infiation model with massive inflaton de- 
caying into another scalar field via quartic coupling [iJl ■ 
Through our simulations, a new and somewhat simpler 
picture of the late stages of preheating emerges. After 
initial transient when instability develops, bubbles form 
and then break-up, the matter distribution soon arranges 
itself in a clumpy state which persists with little changes 
for a long time. One-point probability distribution func- 
tion of total energy density in this state appears to be 
universally lognormal among various two-field preheat- 
ing models. It is tempting to attribute this to relativistic 
turbulence [s^, [13, [4l| . We also see some evidence that 
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the structure formed during preheating continues to grow 
in size on a much longer timescale. Somehow, this picture 
reminds one of large scale structure formation [43I. I4I . I4H . 

This paper is organized in the following way: In Sec- 
tion ini we introduce scalar field models of preheating, 
derive equations of motion, and discuss physical approxi- 
mations we use. Section Hill describes the detailed imple- 
mentation of numerical evolution algorithm. Initial con- 
ditions including quantum fluctuations of the fields pro- 
duced during inflation are discussed in Section llVl with 
particular attention paid to implementation of Gaussian 
random field generator. We briefly recount the theory 
of preheating via broad parametric resonance at the end 
of chaotic inflation in Section |Vl and present our simu- 
lations in Section IVII We conclude by summarizing our 
main results in Section lylll 



II. EQUATIONS OF MOTION 

As our baseline model of reheating we will take a sys- 
tem of N scalar flelds {0*}, minimally coupled to grav- 
ity and interacting through some (non-linear) potential 
V(0*) as described by the action 
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(1) 

The above action can be modified to describe more com- 
plicated geometrical quantities (like a vector field, for ex- 
ample) instead of a real scalar multiplet {0*}. Ultimately 
what we care about is only field equations of motion and 
their gravitational effects, not actual gauge symmetries. 
Variation of the action ([T]) with respect to the metric 
gives Einstein equation with a stress-energy tensor 



T, 



9^u. (2) 



while variation with respect to the field 
of motion for the field 



dV 

w 



gives equation 



(3) 



Although in principle it is possible to solve the com- 
plete system of Einstein and scalar field equations nu- 
merically (see [H, 113, El, li^ for formulation and some 
approaches), in practice, it is not such an easy thing to 
do. Einstein equation solvers in 3 -I- 1 dimensions are 
complex to implement, very expensive to run, and, de- 
spite marked improvement in the recent years (49l |. still 
might have issues with numerical stability. 

Fortunately for us, we do not have to solve the full Ein- 
stein equations. Although preheating is a rather violent 
process, and stress-energy tensor becomes very inhomo- 
geneous due to non-linear field dynamics, the smallness of 
gravitational coupling constant assures that gravitational 
backreaction of these inhomogeneities is rather small (at 



10"'^ level in the simulations presented in this paper). 
Thus we are going to treat the scalar field evolution as 
if it was happening in a homogeneous flat Friedman- 
Robertson- Walker spacetime 



-di^ + a^{t)dx^ 



(4) 



and calculate the inhomogeneous gravitational field due 
to matter distribution with stress-energy tensor Q us- 
ing linear perturbation theory |42| . We will ignore back- 
reaction of metric perturbations on the scalar field evo- 
lution. 

With these simplifying assumptions, the problem be- 
comes much more tractable: we just need to solve a sys- 
tem of coupled scalar field equations of motion ([S]) , which 
in spacetime with metric (|4]) become 



(5) 



Here and later, dot denotes the derivative with respect 
to time t, and spatial differential operators (gradient V 
and Laplacian A) are taken with respect to comoving 
coordinates and three-dimensional flat metric. The ex- 
pansion rate H = a/a and acceleration a are determined 
by averaged Einstein equations 



(p + 3p). 



(6) 



where energy density p and isotropic pressure p are com- 
ponents of the stress energy tensor ([2]) given by 
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6a2 



vm, (8) 



and averages are taken over the whole simulation volume. 

To solve equation ([S]) numerically, one needs to know 
the Hubble parameter value H . Rather than attempting 
to resolve the constraint equation ([6]) at every time step 
(which would result in an implicit evolution scheme) , it is 
faster and more convenient to use the evolution equation 

H = -H^ -^{p + ?,p)^-\y+p) (9) 

to evaluate its value in the future. This is what is done 
in LATTICEEASY. However, we have one more trick up 
our sleeves: a disproportionately huge gain in numerical 
accuracy can be realized by evolving the Hubble length 
L = 1/H instead of the Hubble parameter H by using 



(10) 



For constant equation of state p = wp, Hubble parame- 
ter evolves as if oc 1/t, while Hubble length evolves as 
L (X t. The latter variable has vanishing second (and 



higher) derivatives and correspondingly smaller trunca- 
tion error when discretized to second order. As an added 
bonus, the spatial gradients (which are the single most 
expensive thing to calculate) cancel out when taken in 
p + Sp combination, and do not enter evolution equation 
in the form (fTO)) . 

Once the field equations of motion ^ arc solved, we 
can evaluate all the components of stress-energy tensor 
([2]), in particular energy density ^ and pressure (O, as 
well as calculate linear metric perturbations they create 
in a homogeneous spacetime ([4]). In this paper, we will 
focus on scalar perturbations, behaviour of which dur- 
ing preheating have not been widely studied yet. In the 
longitudinal gauge, perturbed metric can be written as 



(1 + 2$)<it2 + a^(^i _ 2*)rfx2 



(11) 



The two scalar gravitational potentials <& and are equal 
in the absence of anisotropic stress. This is in general 
not the case for the scalar field models, and can be ex- 
pected to hold only approximately and in the average 
sense. Both gravitational potentials <f> and 'J are non- 
dynamical, being solutions of the Poisson-like equations. 
In particular, the equation for 'I' is 



^ T P 



(12) 



Strictly speaking, what should stand in the right hand 
side is not the density p, but the gauge-invariant density 
variable /9,„, which has some extra terms in it |42l |. At 
this stage of the code development, I will not make that 
distinction and simply ignore the missing terms (which 
usually works well on sub- horizon scales), along with 
anisotropic stress contribution to potential $. Full sup- 
port of gauge-invariant perturbations including tensor 
and vector modes is planned for the future code release. 



III. PDE SOLVER IMPLEMENTATION 

Scalar field equations of motion ([5]) are coupled non- 
linear partial differential equations, and have to be solved 
numerically. Fortunately, all the non-linearity comes 
from the potential term only; the differential operator it- 
self is simple, homogeneous and hyperbolic, which makes 
the numerical solution quite straightforward. For the 
solver, we adopt a second-order accurate finite difference 
scheme based on leapfrog algorithm. 

The scalar field values 4>^ are discretized on three- 
dimensional cubic nxnxn grid in comoving coordinates 
with uniform spacing dx and periodic boundary condi- 
tions. Since evolution equations ([5]) are second order, 
values of the fields on two consecutive time slices are re- 
quired to advance to the next one. We will denote the 
previous, current, and next time slices by indices dn, hr, 
and up correspondingly. The time derivatives of a quan- 
tity X arc discretized to second order as 



X 



2dt 



^ _ - 2Xhi- -I- Xdn „x 




FIG. 1: Three-dimensional spatial discretization stencil. 
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TABLE I: Summary of spatial discretization schemes. 

Discretization of spatial differential operators 
D[X] G[X] 



AX 



{dxf 



(vxy 



{dxf 



(14) 



allows more freedom. Direct generalization of the sec- 
ond derivative discretization in equation ([T3|) to three 
spatial dimensions leads to the often-used second-order 
accurate expression for Laplacian using six nearest neigh- 
bours of a point. However, this is not the only (or the 
best) choice. Truncation error for this scheme depends 
on direction, introducing anisotropic artifacts in the field 
evolution at short length scales. One can mitigate this 
unpleasant fact by increasing resolution, but smarter dis- 
cretization scheme is a better solution. By using all 26 
neighbours in a 3 x 3 x 3 cube around a point, one can 
derive a family of discretizations of Laplacian operator 
which is second-order accurate and fourth-order isotropic 
[50l |. For discretization to be isotropic, the coefRcients in 
a linear combination I^fX] approximating the Laplacian 
operator should only depend on the distance from the 
central point 



x+l y+1 z+1 
x—1 y—1 z—1 



(15) 
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as illustrated in Figure [T] The values of coefficients 
for possible discretizations of the Laplacian operator are 
summarized in Table [H along with their computational 
cost and stability properties. Isotropic discretizations A 
and B offer reduced computational cost, while isotropic 
discretization C has the best accuracy and stability. As 
multiplications are cheap and additions are essentially 
free on modern CPUs, there is no reason not to use the 
best discretization scheme available. Thus, DEFROST 
is configured to use the isotropic discretization C by de- 
fault. That can be changed by uncommenting the other 
coefficient definitions in the source code, although profil- 
ing shows little gain in doing so for large grids (for which 
performance of DEFROST solver is apparently memory- 
bandwidth dominated). 

To calculate the energy density ([7]), we need to dis- 
cretizc the square of spatial gradients as well. It is very 
important for discretized energy to be conserved by dis- 
cretized equations of motion. Otherwise, it will leak off 
the grid in the course of a long simulation, affecting over- 
all accuracy or even giving incorrect results (for equa- 
tion of state, for example). Simply squaring the first 
spatial derivative discretized like in equation (|13|) . al- 
though second-order accurate, is not conservative. Us- 
ing discretized action approach [5l[, one can show that 
the conservative second-order accurate and fourth-order 
isotropic discretization of the gradient-square operator is 



fields are accumulated to evaluate 



x+l y+1 2+1 



G[X] 



(16) 
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where coefficients are the same as in Laplacian (|15p. 
Evaluating this expression requires 30 multiplies per 
point for discretization C, which is significantly more ex- 
pensive than computing the Laplacian. However, this 
does not place much of a burden on the total runtime. 
As we mentioned before, we eliminated gradient-square 
terms from evolution equations, so they only need to be 
calculated for output, which happens much less often. 

Putting everything together, we end up with the fol- 
lowing evolution scheme: the discretized field values are 
advanced to the next time slice using 



2 + D/{aaf - Mfdi^ 
1 + I Hdt 



^Hdt 



^hr " 



1 



■Hdt 



^dn, (17) 



where D is the discretized Laplace operator p5p and a = 
dx/dt. The time step dt has to be small enough both to 
satisfy Courant stability condition (a > const listed in 
TableU]) and to resolve the period of the fastest oscillating 
field. All the coefficients in expression (fT7)) except for the 
effective mass of the i-th field 



M? EE — 



1 dV 



(18) 



are constant on the grid and the same for all fields, and 
thus can be pre-computed outside the evolution loop. In- 
side the same loop, kinetic and potential energy of the 



(p + 3p)hr = -3 H 



x,y,z 



2dt^ 



which is used to advance the Hubble length 



(19) 
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2dt. (20) 



To avoid weak numerical instability associated with even- 
odd slice decoupling, we use Lhr = (^up + ^dn)/2 instead 
of Ljir in the above equation, which then can be solved 
for either as an exact quadratic or iteratively (we use 
the latter in the code). 

Equations (fT7|) and (j20|l provide a complete recipe how 
to advance the field variables to the next time step. Once 
in a while (at user's request) we would also like to calcu- 
late, analyze, and output for visualization some auxiliary 
variables like energy density p and gravitational poten- 
tial Energy density ([7|) and pressure Q are easy to 
calculate once we know the field values and their gra- 
dients (|16p . Finding the gravitational potential is a little 
less trivial, as we need to solve the Laplace equation (fT2|) 
with periodic boundary conditions. The fastest way to 
do it is to use a fast Fourier transform (FFT). Applying 
FFT to the discretized Laplacian operator p5|) . we end 
up with an algebraic equation for ^' in Fourier domain 



2 P(cos — A:^:, cos ^ifc^, cos —k, 



(21) 



where fc € [0, n) is the integer Fourier mode wave- number 
and polynomial P{i,j, k) follows from discretization ([T| 



P{i,j,k) ^ c^,+ii{i+j+k)+b2{i3+ik+ik)+c^ijk, (22) 

with Cfe = 2'^Cfc. Here Ck are once again the coefficients 
of discretization (fTSl) with values listed in Table |T1 

DEFROST implements the evolution scheme (fT7|) in 
Fortran 90, fully taking advantage of capabilities of mod- 
ern hardware and compilers by using both automatic vec- 
torization (over field variables) and automatic paralleliza- 
tion (of the evolution loop). During code development, 
profiling showed that physical memory layout has an un- 
expectedly large impact on performance, which became 
apparent as the solver loop got optimized. Let us briefly 
discuss the storage model which was adopted after some 
investigation. The fields are sampled and stored in 
a large multi-dimensional array smp (N , 0:p,0:p,0:p,3). 
The first index (minor in Fortran index ordering) enumer- 
ates the TV fields. The next three indices enumerate the 
three dimensions of the spatial grid, padded to (p -I- 1)'^ 
elements for reasons discussed below. The last index (ma- 
jor in Fortran index ordering) enumerates the three time 
slices used in the evolution scheme. Rather than copying 
large amounts of data around, the allocation of indices 
for dn, hr, and up slices cycles every time step. 
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^ periodic boundary conditions ^ 





FIG. 2: Memory layout of the spatial grid and field samples. 



The layout of a single time slice is shown in Figure [2] 
(with one spatial dimension suppressed for clarity). The 
n X n X n simulation volume is surrounded by a sin- 
gle cell wide boundary layer, introduced to implement 
the boundary conditions without conditional logic in- 
side the evolution loop. It also allows for an easy tran- 
sition to parallel cluster implementation using MPI. as 
required buffers are already in place. Somewhat counter- 
intuitively, padding the grid by a few extra empty cells 
can significantly improve the runtime for large grids. The 
root cause for this phenomenon probably lies in some in- 
teraction between memory access pattern and hardware 
memory cache algorithm. To get the best FFT perfor- 
mance, the grid size n is usually taken to be a power of 
two. While evaluating p7|) on an un-padded array, mem- 
ory would be accessed with a power of two stride, which 
could conceivably interfere with caching and prefetching 
done by the memory subsystem. Whatever the cause, 
padding the array so that its size p + I is prime (or a 
product of a few large primes) can improve the runtime, 
so the user is advised to experiment. 

Finally, a few words should be said about statistical 
estimators used to analyze the simulation data. The 
power spectral density (PSD) estimator is a fairly con- 
ventional one implemented using FFT. It employs anti- 
aliasing with fourth-order polynomial kernel when fold- 
ing the spectrum into wave-number bins to reduce sam- 
pling noise. The implementation of probability density 
function (PDF) estimator is less conventional, and does 
not use histogram binning at all. Instead, PDF is de- 
rived from cumulative density function (CDF) which is 
obtained by partially sorting the data cube into n quan- 



tile brackets. Although more expensive and harder to 
parallelize, this approach is more robust and offers uni- 
form sampling noise across the distribution. 



IV. INITIAL CONDITIONS 

At the end of the inflation, most of the energy is still 
stored in the inflaton, and all the fields are homogeneous 
except for small quantum fluctuations. But the presence 
of these quantum fluctuations in the fields is essential to 
trigger the dynamical instability leading to preheating. 

Initial conditions in the homogeneous field components 
depend on the model of inflation and are treated as an ex- 
ternal input by DEFROST. They arc straightforward to 
obtain by following the expanding Friedman-Robertson- 
Walker solution during inflation either analytically (if 
possible), or using any of the available numerical ordinary 
differential equation integrators. As inflationary trajec- 
tory is an attractor [s^ , it is easy to find and no partic- 
ular care is needed on where to start tracing it. The full 
three-dimensional simulation of preheating should take 
over at some time near the end of the inflation, but as 
there are other factors at play (such as limited spatial dy- 
namic range available to simulation), the exact moment 
is best decided on case by case basis. 

To make a concrete example. Figure [3] shows expansion 
history of a typical chaotic inflation model with massive 
inflaton V{(p) = to^0^/2 (discussed in more detail in the 
next section). As the Universe is inflating, its horizon 
size L = 1/H stays relatively constant, while the physi- 
cal wavelength of comoving modes grows with the scale 
factor a. The modes which were originally inside the hori- 
zon expand and leave the horizon during inflation. Even- 
tually, inflation ends and the horizon size starts growing 
faster than the scale factor (for instance, L oc a'^/^ during 
matter domination), at which point the modes begin to 
re-enter the horizon. The moment in time when comov- 
ing modes stop leaving the horizon and begin re-entry 
can be taken as the end of the inflation. This happens 
when 



dlnL _ ^ 
din a ' 
or in terms of a slow roll parameter 



H 



I. 



(23) 



(24) 



For simulations presented in this paper, we start exactly 
at the moment when inflation ends ([M)) . and select the 
comoving box size of ^ = IQ/m to cover all the length 
scales of interest, as illustrated in Figure [3l 

As we already mentioned, it is crucial to include quan- 
tum field fluctuations in the mostly-homogeneous initial 
conditions for the preheating instability to develop. The 
spectra of quantum field fluctuations are determined by 
effective masses of the fields involved. Let us briefly re- 
count the standard derivation 0] to establish notation. 
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FIG. 3: Expansion history of the universe in chaotic inflation. 
Numerical simulation starts at the end of inflation (e = —1). 
The size of comoving simulation box is selected to cover the 
scales corresponding to horizon size (red), inflaton mass (thin 
black line), and the wave modes in unstable band (not shown). 



Canonically normalized massive field with Lagrangian 

1 



(25) 



is quantized in flat Minkowski spacetime by promoting 
the field value and field momentum 



dijp 



(26) 



to quantum operators (p and tt obeying canonical com- 
mutation relation 



[^(x),^(x')]=z<5(x-x'). 



(27) 



In flat spacetime, the field operator can be represented 
using Fourier mode decomposition as 



<^(x^) 



1 



d^k r 



(28) 



(27r)S J V2w 

where the mode creation-annihilation operators obey 

[a^,a+]=,5(k-q), (29) 

which directly follows from canonical commutation rela- 
tion (|27p . The mode frequency lo of the massive field is 
related to its wavevector k by a simple dispersion relation 
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(30) 



Using mode decomposition (|28p and commutation rela- 
tion ([29|l . it is easy to show that the two-point correlation 
functions of the field value and momentum are 



(^(x)^(x')) - 
(7r(x)^(x')) - 



1 



(27r)3 
1 

(27r)3 



fk 

f_k 
~2 



,ik(x- 



UJ e 



ik(x- 



(31) 
(32) 



The spectrum of the field tp fluctuations is simply l/(2a;). 

One can repeat this procedure in the background of 
expanding homogeneous universe. In general, the time 
dependence of the modes will be different, but it turns out 
that quantum fluctuations of massive flelds in de Sitter 
spacetime can be approximated quite well with above 
expressions if one simply replaces the field mass m with 
an effective mass 



9 , 
4^- 



(33) 



We will use this approximation in DEFROST. In addi- 
tion, we will follow the established practice of treating 
quantum operators as Gaussian random variables. This 
is an assumption, but not totally unjustified one. Quan- 
tum modes essentially behave classically after they leave 
the horizon [sst . Even the sub- horizon modes (which are 
initially quantum) get large occupation numbers once the 
prehea ting instability kicks in, and can be treated clas- 
sically [ll|. Thus, we will initialize the field fluctuations 
as a Gaussian random field 



(^(x, t) 



1 



(2^)i 



(fk 



bk cos cut 4- Ck sin ujt 



(34) 



where the complex random operators bk and Ck obey 



(35) 



{hbl,) = {ckcl,)=6ik-k') 



to reproduce the two-point correlation functions 

A lot of effort has gone into making the realization 
of random field initial conditions in DEFROST as sta- 
tistically accurate as possible. The straightforward and 
often used way to generate a Gaussian random field on 
a discrete grid is to directly discretize equation ([34|) in 
Fourier space, assigning A:-th mode a random Gaussian 
number with amplitude l/^/2uj. Although simple, this 
procedure is spoiled by the finite grid size effects, and 
does not reproduce correct two-point correlation func- 
tions in the real space [1^ . One ends up with a substan- 
tial lack of power on the scales comparable with the box 
size, which is not surprising if one considers that only a 
few long-wavelength modes "fit" into the box, and naive 
discretization ignores all the power in the infrared part 
of the spectrum which should have been properly aliased 
into those few low-fc modes. 

A wealth of the literature is dedicated to the subject of 
generating Gaussian random fields of a given spectrum, 
particularly in the context of the iV-body simulations of 
the large scale structure [H, M, M, M, M, M,M,^- 
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Wc draw on that experience, and in DEFROST for ran- 
dom field generator we adopt a method described in 
[H, [5§|. Gaussian random field ([M]) with l/(2tj) spec- 
trum is realized by convolving white noise with a spher- 
ically symmetric kernel function 



1 



(27r) 
1 



ikx 



2lu 
k^dk 



(36) 



sin kr 



to: 



kr 



The kernel function ^(r) can be evaluated analytically in 
terms of Bessel functions 



Ki{mr) + 2mrK3{mr) , (37) 



and has a power law ultraviolet divergence ^(r) oc r~2 in 
the limit r — > 0. For the purposes of discretization on a 
finite grid, we have to regularize this divergence, which 
we do by introducing a Gaussian cut-off at some scale q 
below the Nyquist frequency 



1 



k^dk 



sin kr 



(fc2 



2 ^ 

cS) 



kr 



exp 



fc2 



(38) 



The regularized kernel does not have a nice expression 
in terms of elementary functions, but is easy to evalu- 
ate numerically, of course. We use a one-dimensional 
discrete sine transform (DST) on a substantially larger 
grid to calculate it (simply because DST is already pro- 
vided by FFTW libraries we use), but other methods like 
quadrature integrators could be used as well. Once the 
spherically-symmetric kernel ^(r) is evaluated, it is sam- 
pled on the the three-dimensional grid in real space, and 
the random field is initialized convolution 



^(x,0) - 



1 



d^kd^x' bk^{x 



(39) 



(27r)3 

The convolution is implemented using discrete FFT as 

Sfce(x')e*( 

k X 



(40) 



where n is the grid size, dk ~ 27r/£ is the spacing of dis- 
crete wavemodes, and B^. is a complex Gaussian random 
number generated using Box-MuUer transformation 



Bk = v/-21n[/ie 



2iTiU2 



(41) 



from two real random numbers Ui and U2 uniformly dis- 
tributed on a unit interval. The implementation of initial 
velocity generator is entirely analogous (and handled by 
the same procedure), and is not worth repeating here. 

Finally, I have to point out that as of current writing, 
there is a bug in the random number generator imple- 
mentation in LATTICEEASY. The formula (glD is mod- 
ified there to produce two complex numbers using only 
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FIG. 4: Stability bands of Mathieu equation. 



three uniformly distributed real numbers. This results 
in correlated random numbers with quite non-Gaussian 
distribution. Fortunately, the results reported so far do 
not seem to be too much affected by this problem, but 
as non-gaussianity studies are becoming more prominent 
in the modern cosmology literature, some care must be 
taken in proper implementation of random number gen- 
erators. 



V. CHAOTIC INFLATION AND BROAD 
PARAMETRIC RESONANCE 

So far, the discussion was rather general, as DEFROST 
is designed to be easily adoptable to study arbitrary 
models of preheating. This Section describes preheating 
model we selected for first simulations with DEFROST: 
chaotic inflation ending via broad parametric resonance 
0. The development of linear instability in this model 
can be largely understood analytically 0|, and its non- 
linea r dy namics has been widely studied numerically as 
well |12| . [1^ . For all its simplicity, this model has very 
rich dynamics, and still holds surprises. Our very first 
simulations uncover new aspects of the evolution dynam- 
ics in this model, which are reported in the next Section. 

In a minimal form, the model consists of two scalar 
fields: the massive inflaton </> and the massless decay 
product ip interacting via potential 



(42) 



During inflation, the value of the inflaton (f> is large, the 
field is over-damped by the large Hubble friction, and 
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slowly rolls down its potential. As it reaches the value 
of around one in Planck units, the damping dips below 
critical, and homogeneous inflaton starts oscillating with 
decreasing amplitude 



(j>{t) Ri <^{t)smmt, <^{t) 



2 2 

3 mt 



(43) 



Decay field ip is coupled to inflaton, and feels its oscilla- 
tions through modulation of the effective mass; equation 
of motion for the Fourier mode ^pk with wavector k is 

i'k + SHijjk + + 5'*' sin^ m?j i'k = 0. (44) 

If the coupling g is large enough, periodic modulation 
of the field mass leads to strong instability via paramet- 
ric resonance. This can be understood analytically by 
applying general theory of differential equations with pe- 
riodic coefficients Q. If one ignores the expansion and 
the Hubble drag term in equation (jH]), evolution for the 
Fourier mode ipk is given by Mathieu equation 



(A-2gcos27;)V'fc =0, 



(45) 



where we have introduced dimensionless parameters 



A = 2q 



Am? 



(46) 



and time variable 77 = mt to bring the equation into 
canonical form. According to Flouqet's theorem, a gen- 
eral solution of Mathieu equation (|45|) is of the form 
e>^^P(r]), where P{ri) is a periodic function with period tt. 
Floquet exponent /i depends on parameters A and q, and 
there is an elegant way to calculate its value [g^, which 
wc (somewhat reluctantly) will omit here, and just quote 
the final result. Figure H] shows the dependence of Re/i 
on parameters A and g as a density plot. For certain 
parameter values Floquet exponent ^ has positive real 
part, leading to exponential instability of the solution; 
these unstable bands are marked on Figure HI The value 
of A for our problem (I46p is restricted to lie above the red 
line A = 2q in Figure [4l which corresponds to the homo- 
geneous mode with k ~ 0. For sufficiently large coupling 
q ^ 1, large portion of available phase space volume is 
unstable, leading to fast development of instability. This 
regime is known as broad parametric resonance. The 
instability grows on a timescale comparable to 1/m (as 
Floquet exponent values are around Kcfi < 1/3), and 
manifests itself after a few dozen of oscillations of the 
inflaton, which is very fast in cosmological terms. In the 
next Section, we describe the non-linear field evolution 
after this instability develops. 



VI. NUMERICAL RESULTS 

Before we present our simulation results, a few words 
should be said about the units used throughout the code. 



As it is clear from the action ^ , we prefer to work with 
dimensionless scalar fields, rather than canonically nor- 
malized ones. In this convention, one can simply think 
of the values of the scalar fields as measured in units 
of Planck mass mp\. Note that we use reduced Planck 
mass 77ipi = (STrCj^i/^ rather than Afpi = G^^/^. The 
only other scale in the model (jH]) is the inflaton mass 
m. The coupling constant g (which has dimension of 
mass with our scalar field normalization), the frequency 
and wavevectors of the unstable modes, and everything 
else can be referred to it. It is convenient to scale all 
the quantities (except for field values) by the appropri- 
ate power of m, making them dimensionless and suitable 
for numerical analysis. Thus we set m = 1 in the code. 

We simulate the preheating model (|42|) with inflaton 
mass m = 5 • 10~^mpi, the value of coupling constant 
g ~ 100 m in a comoving box of size £ = lO/m (the values 
chosen correspond to the ones used in Ref. [Ijl)- The grid 
size is taken to be 256'^, while the time step dt — 2~^^ /m 
has to be reduced substantially below the Courant limit 
to resolve oscillations of the field tj} (which is initially 
hundred times heavier than inflaton 4>). The simulation 
is started at the end of the inflation when the value 
of the inflaton is « 1.009343 and the Hubble constant 
is H K, 0.50467™. The simulation box is initially about 
five Hubble lengths across. We let the code run until 
t = 256/to, which corresponds to 2^^ time steps. 

To get things going, we first reproduce the previously 
reported results [l2[ on expansion history of the preheat- 
ing model ([42]) . The left panel of Figure[5]shows evolution 
of horizon size during expansion, with a'^/^ growth corre- 
sponding to matter-dominated expansion scaled out. The 
expansion history has a sharp break as instability devel- 
ops, and energy gets deposited into relativistic inhomo- 
geneous modes from a homogeneous oscillating inflaton 
(which behaves as a pressureless dust). This transition 
can be seen in terms of an effective equation of state pa- 
rameter w = (p) / (p) , the value of which is plotted in 
the right panel of Figure [S] It undergoes large ampli- 
tude oscillations (shown by thin pale line in Figure [5]), 
but when averaged over a few periods (with a Kaiser 
window function), the underlying behavior is uncovered. 
The averaged equation of state (shown by thick red line) 
switches over from dust-like equation of state tjj = to a 
value slightly less than a quarter corresponding to 
a fairly relativistic fluid. (If evolved further, the residual 
homogeneous component in the inflaton will eventually 
come to dominate the evolution again, slowly lowering 
equation of state toward w = in the process.) 

While we recover the results obtained by simulations 
with LATTICEEASY, the accuracy of the integrator 
used in DEFROST is significantly higher. The perfor- 
mance of integration scheme for expansion factor (j20p is 
illustrated in Figure [SI which shows residual curvature 



K 



3 



(47) 



which should be zero in the flat model we are evolving. 
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FIG. 5: Expansion history (left) and average 

As you can see, the constraint equation is satisfied to 
10^^ level (which is exceptionally good for a second or- 
der scheme with 10""^ timcstep), and the error does not 
accumulate with time. In fact, this error is mostly due 
to the fact that we neglected second order corrections to 
density from initial field fluctuations. 
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FIG. 6; Residual curvature K testing constraint violation. 
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equation of state (right) during preheating. 

Having made sure that our code reproduces previously 
reported results, and after numerous checks of code in- 
tegrity and accuracy, we move on to investigation of the 
field dynamics during preheating. 

Evolution of field distributions and spectra for infla- 
ton 0, decay field ip, total density p and gravitational 
potential ^' are presented in Figure [71 Left panel shows 
evolution of the median value (thick red line), along with 
68% and 95% percentile brackets around it (which would 
correspond to la and 2a contours for a Gaussian distri- 
bution) shown by shaded outlines. The contour with the 
lightest shading spans the extremal values inside the sim- 
ulation box, which serves to illustrate the extent of the 
tails of the distribution, although the exact percentile it 
corresponds to depends on the spatial resolution of the 
simulation. Dilution due to expansion has been scaled 
out to highlight the relative change of the distributions 
as evolution proceeds. 

The evolution of the distribution of the decay field val- 
ues (second row of Figure [7|) clearly shows the onset of 
instability a little after t = lOO/m, rapid spreading of 
the distribution due to exponential amplification of seed 
inhomogeneitieS; and self-limiting of the growth by non- 
linear interactions when the scaled value of the decay field 
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FIG. 7: Evolution of field value distributions (left) and spectra (right) of inflaton (j) (top row), decay field t/j (second row), total 
density p (third row), and gravitational potential 'i' (bottom row). Onset of instability and characteristic size of the structure 
formed is clearly visible on the plots. 
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becomes of order unity. As the decay field perturbation 
grows and becomes non-linear, it is drawing the energy 
from the zero mode of the inflaton (top row of Figure [7]), 
reducing the amplitude of its oscillations, and eventually 
forcing the inflaton to become strongly inhomogeneous 
as well due to non-linear backreaction. We should note 
here that although the amplitude of the coherent infla- 
ton oscillation decays, it does not go away altogether, 
and the left-over homogeneous mode will eventually come 
to dominate the universe expansion [Tsj . as its equation 
of state is effectively that of pressureless dust, and it 
dilutes slower than inhomogeneous components, which 
have equation of state closer to relativistic one. 

Of particular interest to us is the distribution of the 
total energy density, shown in the third row of Figure [T] 
It is clearly very inhomogeneous, with peak densities eas- 
ily exceeding ten times the average. After a brief tran- 
sient, it quickly settles to a nearly stationary distribution, 
which appears to be highly non-Gaussian (and is in fact 
plotted on a logarithmic scale) . We will come back to this 
point after we inspect the spatial picture of the energy 
density distribution inside the simulation box. 

The bottom row of Figure [7] shows the evolution of 
the gravitational potential. Despite total energy den- 
sity being highly inhomogeneous and having huge over- 
densities, the gravitational potential it produces is rather 
small (at 10"'^ level), and is further diluted away by the 
expansion with near-relativistic equation of state. The 
maximal potential well depth inside the simulation box 
is 2 • 10"'^, which is far too small to form any primor- 
dial black holes. This result is in line with observations 
of [23, [H], although now we have a more transparent 
diagnostic of black hole formation as we calculate the 
gravitational potential directly. 

The right panel of Figure [7] shows evolution of the 
field spectra as the density plot in terms of time t and 
wavenumber k, with dilution of the fields due to expan- 
sion and overall power-law dependence on k scaled out. 
The spectra of total energy density pk and gravitational 
potential ^fc show a very clear simultaneous peak soon 
after instability develops, which is sharply localized in 
both time and scale. Subsequent evolution of the two 
differs, however. The peak power in energy density is 
evolving toward higher k and smaller scales, while the 
peak power in potential is evolving to lower k, so the 
structure of potential wells is growing in spatial extent. 
The understanding of what's going on will become more 
clear when we look at evolution in real space, which is 
what we are going to discuss next. 

Figure [5] shows three-dimensional volume renders of 
the contents of the simulation box. As the fields and 
-0 oscillate rapidly in a standing wave pattern, quickly 
losing coherent phasing, the view of their values is messy 
and not very enlightening, and we will omit it here. Much 
more interesting is the picture of what is happening to 
the total energy density, which is an adiabatic invariant 
for the oscillating fields. The top row of Figure [5] shows 
density distribution inside the simulation box soon af- 



ter onset of instability (at t = 124/m, left) and at the 
end of the simulation (at t = 256/m, right). Density 
field is shaded using logarithmic color map with linear 
transparency ramp applied, so that only the peaks of the 
density distribution are visible. 

Immediately after the onset of instability, the density 
distribution in Figure [5] (top left) looks like smoke is fill- 
ing the box. What you are seeing is actually the over- 
dense bubble walls forming a three-dimensional foam-like 
structure that fills the box. Its origin is easy to un- 
derstand if one thinks about how seed inhomogeneities 
are amplified by instability. Broad parametric resonance 
amplifies wavemodes in a certain band, effectively serv- 
ing as a low-pass filter (with a kernel that can be ap- 
proximated analytically) and sets the characteristic size 
of the structure which grows out of the seed inhomo- 
geneities. Original fluctuations are a Gaussian random 
field, which already has the structure of peaks, ridges, 
and valleys imprinted into it. The skeleton of this struc- 
ture is essentially preserved unchanged as the growth 
of inhomogeneities due to instability increases density 
contrast. Once the density contrast becomes of order 
unity, non-linear evolution takes over. This will happen 
to under-dense regions first, with repulsive interaction 
term g^cfP'ilP helping to evacuate the bubble interiors, 
and pushing the matter density into the bubble walls, 
thus forming the structure you see in Figure [5] (top left). 

The evolution does not stop at forming bubbles, how- 
ever. The repulsive interaction soon breaks the extended 
bubble walls into smaller localized blobs, moving more or 
less freely inside the simulation box (the animation of this 
process is available ' onlineP . The final state is depicted in 
Figured] (top right), and persists with little change for 
a long, long time. This state seems quite distinct from 
thermal equilibrium, yet it is still long-lived and statis- 
tically simple in a certain way, which is quite surprising. 
Even more surprisingly, the distribution of values of total 
energy density p in units of quickly becomes statis- 
tically stationary and after a brief transient tends to a 
distribution with a probability density function shown in 
the lower right of Figure [H It can be fitted with exceed- 
ingly high accuracy by a lognormal distribution 
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(48) 



with one free parameter (tj = 0.6584 or /x = —0.2197), as 
the mean p = exp(^ -I- ct^/2) is unit-normalized by virtue 
of us scaling out the expansion. The corresponding me- 
dian is e'' = 0.8028. With statistical errors of PDF esti- 
mator being what they are. the apparent lognormality of 
a density distribution is undoubtedly not a mere coinci- 
dence, but must have a explanation rooted in scalar field 
dynamics. Moreover, further simulations of preheating 
models with different couplings and inflaton potentials 
seem to suggest that lognormal distribution of density is 
a universal feature of two-fleld preheating models. This 
observation presents a very interesting theoretical puzzle, 
which will be explored in detail elsewhere [631] . 
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FIG. 8: Volume rendering of total density p (top row), gravitational potential 'I' (second raw), and PDF of density values p 
inside the simulation box soon after onset of instability (at t = 124/m, left) and during subsequent evolution (at t = 256/m, 
right). Animations for density, and potential evolution are available at http: //www, sf u. ca/physics/cosmology/def rost/, 
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FIG. 9: Evolution of characteristic structure size in total energy density (left) and gravitational potential (right). 



Although one-point distribution of total energy density 
p quickly becomes stationary as noted above, other quan- 
tities continue evolving on much longer time scales. The 
blobs continue to fragment, and their characteristic size 
slowly decreases with time. While this is obvious from 
visualizations, it can be further quantified by introduc- 
ing the (physical) correlation length of the total energy 
density configuration 



{jap? 



(49) 



The evolution of the comoving correlation length £p/a 
is plotted in the left panel of Figure [51 Initially it is 
very large (> 100/m), as the density field is nearly ho- 
mogeneous. As instability develops and structure forms, 
it abruptly drops to about 10~^/to, and then continues 
to decrease, but much more slowly. Although the graph 
clearly shows evolutionary trend in density correlation 
length £p/a, actual numbers should be taken with a grain 
of salt, as the density field does eventually become frag- 
mented on a scale close to spatial grid resolution (which 
for 256'^ grid we use would be reached around t ~ lO'^/m). 

While it is known that thermalization after preheating 
might take a long time Q , and there might be an interme- 
diate scaling regime in the evolution of the fields [13, EH , 
the view of the process from the real space presented here 
is strikingly simple. One would think the field distribu- 
tions will thermalize eventually, presumably forming a 
homogeneous fluid-like state with Maxwellian distribu- 
tion of particle velocities. Thermalization does not hap- 
pen in our simulations, which lack necessary quantum 
effects. Exactly how and when it does happen is an ex- 
tremely interesting question, and the one which requires 
further study. 

Finally, let us discuss gravitational potential which 
is sourced by the evolving energy density distribution p, 
and is shown in the middle row of Figure [51 To make the 
structure more visible, we have opted for a density plot 
on a three-slice through the simulation box rather than 



a volume rendering. The color map shows positive po- 
tential values (corresponding to under-dense regions) as 
shades of red, and negative potential values (correspond- 
ing to over-dense regions) as shades of blue, blending into 
white for zero potential value. 

Immediately after the onset of instability, the gravita- 
tional potential in Figure [8] (middle left) clearly traces 
the foam-like structure of matter distribution. The iso- 
lated potential peaks (red) in the interior of the bubbles 
are separated by extended potential valleys (blue) cre- 
ated by over-dense bubble walls. The gravitational po- 
tential configuration is asymmetric between positive and 
negative values, and is clearly non-Gaussian. Subsequent 
evolution of the gravitational potential is rather interest- 
ing. As bubble walls break into smaller and smaller blobs, 
the structure of the gravitational potential does not fol- 
low suit. Instead, it begins to grow in spatial extent (the 
animation is available online ) . By the end of the simula- 
tion in Figure[8l (middle right), the size of the structure in 
the gravitational potential spans almost the whole box. 
The growth of structure can be quantified by introducing 
correlation length £qr for gravitational potential the way 
we did for energy density 
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Evolution of the comoving correlation length i^,/ais plot- 
ted in the right panel of Figure [51 with pale line tracing 
its instantaneous value, and thick red line giving a run- 
ning average (with Kaiser window) over few oscillations. 
While comoving correlation length grows overall, evolu- 
tion shown in Figure [51 (right) still represents an initial 
transient, and the long-time asymptotic behaviour is not 
reached in the simulation reported here. Further investi- 
gation shows that correlation length £^^, continues to grow 
faster than comoving size, but not quite as fast as the 
horizon size. Eventually one might even have to worry 
about it outgrowing the finite simulation box size, but 
that would take a very long time, and is not reached by 
our simulations. 
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VII. CONCLUSIONS 

This paper presents a new numerical code I de- 
veloped for simulating preheating of the Universe 
after the end of the inflation, which I call DE- 
FROST. It is small (about 600 lines of Fortran code), 
fast, easy to modify, and is fully instrumented for 
3D visualizations (using [LLNL's Visits for example). 
The source code is available for download online 
at http : //www . sf u . ca/physics/cosmology/def rost/ 
and is distributed under the terms of GNU Public Li- 
cense. While the main design goal of DEFROST has been 
the accuracy of the simulations, performance of the solver 
has also been significantly improved compared to LAT- 
TICEEASY [131, which is the most mature and widely 
used reheating code publicly available today. 

As a result of all the optimizations (and a bit of 
black magic), DEFROST outperforms LATTICEEASY 
by about a factor of four in raw PDE solver speed (for 2 
fields on a 256'^ grid in double precision on a dual Xeon 
5160 machine) while using more accurate (and more ex- 
pensive) discretization. If one takes into account the time 
spent on analysis of the results, the difference is even 
larger, as FFTW libraries used by DEFROST are vastly 
faster than EFT routines shipped with LATTICEEASY 
(especially on multi-processor machines). The speed-up 
offered by DEFROST is so significant that the studies 
done a few years ago on a big parallel cluster [T^, [TJ us- 
ing MPI version of LATTICEEASY [H] can now be car- 
ried out on a single fast workstation. The planned MPI 
version of DEFROST should be able to push the acces- 
sible simulation size over the 1024"^ barrier, provided the 
code scales well. 

The code was tested on a number of chaotic inflation 
models which end via parametric resonance. In this pa- 
per, we report the simulations of the simplest two-field 
preheating model with massive inflaton and quartic cou- 
pling to decay field ([42|) . Wc reproduce the previously 
published numerical results for this model [U [U (and 
the ones for trilinear coupling (isj . which we will not 
discuss here, although our simulation data and results 
for that model are available lonlinel as well) . We further 
investigate the dynamics of the scalar field evolution in 
these preheating models, taking advantage of advanced 
visualization and analysis capabilities DEFROST offers. 
In particular, we study the behaviour of energy density 
distribution and scalar gravitational potential during pre- 
heating, something which has not been looked at closely 
before. Our main science results are summarized by two 
observations, both novel and quite surprising. 

First, the evolving scalar fields quickly end up in a 
simple state, which, although highly inhomogcneous, ap- 
pears to have a certain universality to it. In this state, the 
one-point distribution function of total energy density is 
nearly stationary (apart for the overall dilution due to 
expansion), and is described by a lognormal distribution 
for all two-field parametric resonance preheating models 
we tried so far, namely the ones described by interaction 



potentials 

. !/(</), V^) = lA^-'^igW'", 

. y(0,^) = iA(</>2-HVT. 

This is true even if distributions of field values or other 
correlators might still be evolving, and appears to be a 
very general statement about random scalar fields one 
encounters in preheating. It is tem pting to attribute this 
state to scalar field turbulence [40, [4l| , especially since 
lognormal density distributions are known to occur in 
supersonic isothermal turbulence in hydrodynamics [soj . 
We do not see obvious signs of thermalization, even if 
the simulations are run for a time much longer than the 
dynamical timescale of the problem (the longest done so 
far for massive inflaton is corresponding to five e- 

folds since the end of inflation; this is limited mainly by 
my patience rather than the code stability). 

Second, less general but still amusing, is the observa- 
tion that the small-scale structure in the gravitational 
potential can grow faster than comoving box expands. It 
is not quite clear whether the reason for it happening is 
kinematic or dynamical in nature. As we neglected grav- 
itational interactions in our simulations, the only thing 
that can cause the structure to grow is the interaction be- 
tween scalar fields themselves. In our preheating model 
it is repulsive; yet the structure still grows! Although 
one might suspect that any inhomogeneity in gravita- 
tional potential on sub-horizon scales would probably get 
washed away by subsequent evolution (and is too small to 
form primordial black holes), this effect still might have 
some interesting cosmological consequences. 

All in all, we find that the picture of preheating dy- 
namics is simpler in real space than what it looks like 
in particle representation. The final stage of preheating, 
with growing structure and lognormal density distribu- 
tion, eerily reminds one of large scale structure forma- 
tion in later cosmology (although of course it occurs on 
vastly smaller scales and is driven by completely differ- 
ent physics). Perhaps the analytical methods developed 
for the latter [H, O, \^ could be fruitfully applied to 
preheating as well. This is what we intend to explore 
next. 
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